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ABSTRACT 


Stochastic robustness , a simple technique used to estimate 
the robustness of linear, time-invariant systems, is applied to a 
single-link robot arm control system. Concepts behind 
stochastic stability robustness are extended to systems with 
estimators and to stochastic performance robustness. 
Stochastic performance robustness measures based on classical 
design specifications are introduced, and the relationship 
between stochastic robustness measures and control system 
design parameters is discussed. The application of stochastic 
performance robustness, and the relationship between 
performance objectives and design parameters are demonstrated 
by means of the example. The results prove stochastic 
robustness to be a good overall robustness an' 'is method that 
can relate robustness characteristics to con. stem design 
parameters. 

INTRODUCTION 

Standard linear control system design techniques rely on 
accurate models of the system to be controlled. Because 
models are never perfect, robustness analysis is necessary to 
determine the possibility of instability or inadequate 
performance in the face of uncertainty. Robustness to these 
uncertainties, parametric or unstructured, is normally treated 
deterministically and often without regard to possible physical 
variations in the system. Consequently, overconservative 
control system designs, or designs that are insufficiently robust 
in the face of real world uncertainties are a danger. 

Stochastic robustness, a simple technique to determine the 
robustness of linear, time-invariant systems by Monte Carlo 
methods was introduced in [1] and presented in detail in [2,3]. 
These references described stochastic stability robustness 
analysis and introduced the probability of instability as a scalar 
measure of stability robustness. Confidence intervals for the 
scalar probability of instability were presented, and the 
stochastic root locus, or probability density of the closed-loop 
eigenvalues, has shown to portray robustness properties 
graphically. Because it is a statistical measure of robustness, 
and because it directly uses knowledge of the statistics of the 
physical parameter variations, stochastic robustness is 
inherently intuitive and precise. The physical meaning behind 
the probability of instability is apparent, and overconservative 
or insufficiently robust designs can be avoided. Applications 
of stochastic robustness to analyzing full-state feedback aircraft 
control systems were described in (4J. The results presented 
there illustrated the use of stochastic stability robustness 
techniques in comparing control system designs and in 
including finite-dimensional uncertain dynamics. 

Concepts behind stochastic stability robustness can be 
easily extended to provide insight into control system design 
for performance. Design specifications such as rise time, 
overshoot, settling time, dead time, and steady-state error are 
normally used as indicators of adequate performance and lend 
themselves to the same kind of analysis as described above. 
Concepts of stochastic stability robustness analysis can be 
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applied to these criteria giving probabilistic bounds on 
individual scalar performance criteria. Metrics resulting from 
stability and performance robustness can be related to controller 
parameters, thus providing a foundation for design tradeoffs 
and optimization. Details of these extensions and uses of 
stochastic robustness are described in the following and are 
illustrated by means of an example. 

STOCHASTIC STABILITY ROBUSTNESS 

Stochastic stability robustness of a linear, time-invariant 
(LTI) system was described in [3] and is summarized here. 
Consider a LTI system subject to constant-coefficient control: 

x(r) = F(p)x(r) + G(p)u(0 (1 

y(r) = H(p)x(r) (2 

u(f) = ii r (f) - CH(p)x(0 (3 

x(/), u(0. y(0i and P are state, control, output, and parameter 
vectors of dimension m, q> and r, respectively, accompanied 
by conformable dynamic, control, and output matrices that may 
be arbitrary functions of p. u c (r) is a command input vector, 
and, for simplicity, the (m x n) control gain matrix C is 
assumed to be known without error. The n eigenvalues, X/ = 
Oi + jeo/, i = 1 to n, of the matrix [F(p) - G(p)CII(p)] 
determine closed-loop stability. The control gain matrix C is 
designed using some nominal or “mean” value of the dynamic 
model, denoted F, G, and H, that represents F(p), G(p), H(p) 
evaluated at the nominal parameter vector. The actual system 
has an unknown description, denoted F a, Ga, and IIa that 
depends on the actual (unknown) value of the parameter vector 
p. Environment, variations in the nominal state, system 
failures, parameter estimation errors, wear, and manufacturing 
differences all can contribute to mismatch between the actual 
system and that used to design the controller. The parameter 
vector p is assumed to have a known or estimated probability 
density function, denoted pr(p), that expresses the statistics of 
parametric uncertainly due to the above factors. 

System stability requires that no eigenvalues have positive 
real parts. While the relationship between parameters and 
eigenvalues is complicated, estimating the probability of 

A 

instability (F) of a closed-loop system from repeated eigenvalue 
calculation is a straightforward task. Using Monte Carlo 
evaluation , the closed-loop eigenvalues are evaluated J times 
with each element of py, j = 1 to 7, specified by a random- 
number generator wdiose individual outputs are shaped by 
pr(p). The probability- of" -instability estimate becomes 
increasingly precise as 7 becomes large. Then, 

N(amax £0) 

Instability) = lim ~ “ -y (4 

y _> oo J 

and 

Pr(instability) = F = 1 - Pr(stability) (5 
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/V(.) is the number of cases for which all elements of a , the 
vector of the real parts of the closed-loop eigenvalues, are less 
than or equal to zero, that is, for which a m ax ^ 0, where 
Gmax is the maximum real eigenvalue component in a. For 
less than an infinite number of evaluations, the resulting Monte 

A 

Carlo evaluation is an estimate, denoted P. 

Because P is a binomial variable (i.e., the outcome of 
each Monte Carlo evaluation takes on one of two values: stable 
or unstable) confidence intervals are calculated using the 
binomial test, where lower (L) and upper (U) intervals satisfy 
[5] 

Pr(X < x-1) = ^(7,;) L/(l - L) ; - J = I - j (6 
7=0 

Pr(X < x) = ZAJ, j ) \jJ ( 1 - \J) J ' J = j (7 

7=0 

X is the actual number of unstable cases after J evaluations (X 
a yi . . 

-J P) and (7, x) is the binomial coefficient, - y- Explicit 

approximations of the binomial test [6, 7] avoid an iterative 
solution of Eq. 6 and 7 for (L,U) and are accurate to within 
0. 1 %. Confidence intervals for the unknown, true probability 
of instability are presented in more detail in [3]. 

Stochastic Stability Robustness of Systems with 
Estimators 

Stochastic stability robustness analysis is easily extended 
to systems incorporating dynamic state estimators. Using Fa, 
Ga, and Ha as the actual system matrices and F, G, and II as 
the design system, state and estimator equations are [8] 

x(0 = F a (p)x(0 + G a (p)u(0 + L a (p)w(0 (8 

*(r) = F*(r) + Gu(0 + K[ z(r) - H*(i)| (9 

u(o = u c (0 - c$(o (in 

where $(r) is the estimated state, L is the disturbance input 
matrix, and K is the estimator gain matrix. The measurement is 
taken through the actual output matrix and can include 
measurement noise: 

z(r) = II A (p)x(0 + n(0 (11 

Stability of the full dynamic compensator is determined by the 
combined system. Substituting for the control in Eq. 8 and 9, 
the eigenvalues of the closed-loop system matrix F c i must be 
in the left-half plane for stability: 

[F a -GaC 
hc! "LKH A F-GC-KH 

In terms of the state (x) and error dynamics (x - £), the 
coupling effect of mismatch on the closed-loop system is 
directly apparent (8] 

[ F a -GaC -GaC 1 

Fcl “ L(F-Fa) - (G-G a ) - K(IMIa) F - (G-Ga)C - KIlJ 

(13 

The effect of parametric uncertainty on stability robustness is 
computed by Monte Carlo evaluation of the eigenvalues of Eq. 
10, with F(p), G(p), and H(p) substituted for F A , Ga, and 
Ha- Closed-loop eigenvalue densities portrayed on the 



stochastic root locus show the possible interaction of dynamic 
and estimator states, and the possible robustness degradation 
due to the estimator. Well-known loss of LQ stability margins 
when a state estimator is added [9] can be quantified by the 
probability of instability. 

STOCHASTIC PERFORMANCE ROBUSTNESS 

While stability is an important element of robustness 
analysis, performance robustness analysis is vital to 
determining whether important design specifications are met. 
Stochastic stability robustness is described by a single 
parameter, the probability of instability. Adequate 
performance — initial condition response, response to 
commanded inputs, control authority, and rejection of 
disturbances is difficult to describe by a single scalar. 
However, elements of stochastic stability robustness analysis 
(e.g., Monte Carlo evaluation and use of the tests described 
above for confidence intervals) apply independent of the 
performance criteria chosen. This will be demonstrated in the 
sequel. 

Numerous criteria stemming from classical control 
concepts exist as measures of adequate performance. 
Appealing to these, one can begin a smooth transition from 
stability robustness analysis to performance robustness analysis 
simply by analyzing the degree of stability or instability 
rather than strict stability. As described in [3], one method of 
doing this is to shift the vertical discriminant line from zero to £ 
less than (or greater than) zero (Fig. la). Histograms and 
cumulative distributions for varying degrees of stability are 

readily given by the Monte Carlo estimate of 1 - fpr(o)da, 

- oo 

w'here £ represents a maximum real eigenvalue component, 
and -eo < £ < oo. The histogram is a plot of 

N [(I - A) <P jji ax_gJLL vs £. £ is an increment in £, N\. ] is 

the number of cases whose maximum real eigenvalue 
components lie in the increment, and J is the total number of 
evaluations. The histogram estimates the stability probability 
density function , pr(£), which is obtained in the limit for a 
continuous distribution of £ as A -> 0 and J -> °°. The 
cumulative probability distribution of stability , Pr(£), is 

similarly estimated and presented as — — m j x * vs. £, the 

exact distribution being achieved in the limit as J -> 
Binomial confidence intervals are applicable to each point of the 
cumulative distribution, as there are just two values of interest, 
e.g., ‘‘satisfactory” or “unsatisfactory”. The probability of 
instability is a special case where £ equals zero. 



Figure I Alternate regions for probability of closed-loop eigenvalue 
location to aid in robustness analysis of transient response. 

a) Alternate discriminant 

b) Sector defined along lines of constant damping and natural frequency. 

The measure of robustness resulting from the cumulative 
probability distribution is directly related to classical concepts 
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of rates of decay (growth) of the closed-loop response, time-to- 
half and time-to double: 


_ 0-693 

f>n 


(14 


0.693 

t double - - £ Wn 


(15 


input, disturbance, initial condition, or some combination is 
evaluated J times, and for each evaluation, the trajectory is a 
binomial variable; it either stays within the envelope or violates 
the envelope. Although computing actual time responses is 
more computation-intensive than probabili ty-of-instability 
estimation or estimation of scalar values associated with 
performance robustness, such analysis is well within the 
capability of existing workstations. 


Taking degree-of-stability analysis one step further, rather than 
a vertical discriminant line, one can confine the closed-loop 
roots to sectors in the complex plane bounded by lines of 
constant damping and lines of constant natural frequency (Fig. 
lb). Roots confined to these sectors would be expected to 
display a certain transient response speed. Again, the 
probability of roots lying within a sector is a binomial variable, 
and confidence interval calculations presented above apply. 

While the speed of the transient response depends on the 
closed-loop poles, its magnitude and overall shape depends on 
the coefficients of the characteristic exponential and sinusoidal 
terms. Closed-loop zeros, residues, eigenvectors, and steady- 
state response are all concepts related to the magnitude of the 
response. The distribution of closed-loop zeros, of residues, 
of steady-state response, and of important elements of the 
eigenvectors all can be estimated by Monte Carlo analysis. 
Hence, a qualitative idea of the possible closed-loop responses 
can be obtained without calculating actual time responses. 

Time responses provide the most clear-cut means of 
evaluating performance, but they are the most computation- 
intensive means as well. If actual time-responses are 
computed, stochastic performance robustness can be portrayed 
as a distribution of possible trajectories around a nominal or 
desired trajectory. After defining “envelopes” around the 
nominal trajectory (Fig. 2), the probability of violating the 
envelopes can be computed using Monte Carlo evaluation. The 
envelope chosen around the nominal trajectory encompasses 
scalar performance measures; the trajectories in Fig. 2 are 
examples of bounds defined by minimum and/or maximum 
allowable dead time, delay time, rise time, time to peak 
overshoot, peak overshoot, settling time, and steady-state 
error. While it is simple to conclude that a response violates 
such an envelope, individual responses within the envelope 
may not be acceptable. In such cases, the derivative of a 
response and envelopes around the derivative also can be used 
as performance criteria. 


min max peak overshoot max 



max 

undershoot 

Figure 2. Example of step response bounds formed by scalar performance 
characteristics. 

There is no unique set of criteria defining envelopes that 
bound an acceptable time response; the segmented envelopes in 
Fig. 2 can be smoothed, or other scalars can be used to define 
points on the envelope. However, once an envelope is 
defined, time response distributions can be computed by Monte 
Carlo methods. The closed-loop time response to a command 


STOCHASTIC ROBUSTNESS AS A CONTROL 
DESIGN AID 


Stochastic robustness metrics can and should be related to 


control system design parameters for robust control system 
design. While general “rules of thumb” regarding the design of 
robust control systems are useful, stochastic robustness metrics 
can identify non-obvious robustness “structures” of particular 
applications. Figure 3 illustrates tradeoffs that can exist and be 
uncovered by stochastic robustness analysis. Consider Fig. 
3a, which shows the upper-half plane of a plant that has a 
complex pair of poles and a right-half-plane zero. Hypothetical 
“uncertainty circles” are drawn around possible closed-loop- 
root locations. As gain increases along the root locus, the 
uncertainty is magnified, and stability robustness decreases, 
w ith possible closed-loop root locations in the right-half plane 
at high enough gain. This case illustrates one where the 
decrease in robustness may be monotonic. Figure 3b 
postulates a system with a real pole and a complex pair of poles 
and zeros. In Fig. 3b, the hypothetical root locus nears the jeo 
axis before ending at the zero in the left-half plane. Again, 
uncertainty circles enlarge as gain increases. In this case, it is 
possible that the root distributions cross into the right-half 
plane, yet are entirely in the left-half plane as gain increases 
further. Here, stability robustness (as measured by the 
probability of instability) may have local or global minima (as a 
function of gain). For multivariable systems with many 
parameters, the intrinsic structure of the problem and the 
tradeoff between the spread in closed-loop-root-location 
uncertainty versus the magnitude of the control gains may not 
be immediately evident. Plots of stochastic robustness metrics 
versus scalar control design parameters provide the necessary 
insight. 



Figure 3 Illustration of design insight revealed by stochastic robustness 
analysis. Solid points indicate closed-loop eigenvalues enclosed by 
“uncertainty circles” . 

a) Root locus sketch where stability robustness decreases monotonically 
with increased gain. 

b) Root locus sketch where stability robustness decreases then increases 
with increased gain. 


EXAMPLE: SINGLE LINK ROBOT ARM 

As an example of the concepts discussed above, stochastic 
robustness is applied to the flexible one-link robot described in 
Ref. 10. The model retains the first three flexible modes, and 
the tip of the link is controlled by applying a control torque to 
the hub, or base of the link. Because the structure of this 
model is representative of a general flexible structure, and the 
physical parameters are easily identifiable, it is a good candidate 
for stochastic robustness analysis. 

The dynamic, control effect, and output matrices are given 
by 
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0 10 0 0 0 

0 0 0 0 0 0 

0 0 0 1 0 0 

0 0 ( 1 )^ - 2 ^ 10)1 0 0 


0 0 
0 0 
0 0 

0 0 


F = 


0 0 
0 0 


0 0 
0 0 


0 1 0 0 

«2 - 2 C 2™2 0 0 


Degree of stability is portrayed in Fig. 5 by the histogram 
and cumulative distribution of stability around the origin for a 
bin size of 0.25. (The smaller bin size was chosen for better 
resolution.) For binary parameter variations [3] of the same 
magnitude as the maximum uniform variations, 2^ or 16,384 
deterministic evaluations give a zero probability of instability as 
well. These results indicate good stability robustness in the 
face of reasonably large uncertainties. 
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(17,18 


where .t is the length along the arm, *i(*),are the normal 

modes, *' = L is the length of the arm, and Ij is the total 

inertia of the arm. The measurements taken through II are the 
tip displacement and hub-rate, respectively. The flexibility of 
the open-loop system is apparent in the open-loop eigenvalues, 
which are 0, 0, -0.177 + 1 1 .8 1 j, -0.432 + 2 1 .6 1 j, and -0.968 
+ 48.37j. The system has a readily identifiable 14-element 
parameter vector: 

P = Kl «l C2 W2 C3 €03 *1(0) *2(0) *3(0) L ML) *2 (L) *3 W h \ 

(19 

Details concerning the modeling and parameter identification are 
given in [10]. The linear-quadratic-regulator designed in [10] 
was used for demonstration and review of stochastic stability 
robustness analysis. The quadratic performance index weights 
tip position and tip-rate, and the LQR state-weighting and 
control-weighting matrices are 

Q = 0.01F t II t [o o] i,f + !,t [o 2]h.r = 0.001 

( 20.21 


resulting in linear-quadratic regulator gains 

C = [35.42 13.38 41.24 2.65 59.32 -0.67 135.46 1.581 (22 

and closed-loop eigenvalues -5.41 + 48. 8j, -6.47 + 23. 8j, 
-6.1 + 2.66 j, -7.7 ± 1 1.42j. 

A uniform probability density function was chosen to 
model the parameter uncertainty statistics, with variations 
between +2% of the nominal values for L and If and +25% 
for the remaining parameters. The stochastic root locus for 
50,000 evaluations of the full-state feedback case is given in 
Fig. 4. The nominal eigenvalues are marked, and the possible 
distribution is indicated by the height above the complex plane 
in units of roots/length along the real axis and roots/area in the 
complex plane. The “bin” size in Fig. 4 is 0.9 along the axis 
and 0.9 x 0.9 off of the axis. The probability of instability 

estimate (?) for 50,000 evaluations is zero. Each of the four 
complex eigenvalue pairs appears in Fig. 4 as a “peak”, with 
some distribution around the peak due to parametric 
uncertainty. The peaks can be well-defined (as in the lowest 
frequency complex pair) or broad (as in the highest frequency 
pair), and the nominal eigenvalues are not necessarily at the 
peaks of the distributions. Uncertainties of the magnitude 
chosen cause complex pairs to coalesce into real roots resulting 
in a distribution along the real axis. 


Moving to performance robustness analysis, Fig. 6 shows 
a top, contour-shaded view of a stochastic root-locus with 
sector bounds defined by 4 < 03 n < 65 and £ > 0.1. For 
50,000 evaluations, the probability of having closed-loop 
eigenvalues outside of these bounds is 0.0147, with 95% 
confidence intervals (L,U) = (0.0136,0.0158). While the 
shape of the time response depends on closed-loop zeros, a 
minimum speed of response can be guaranteed by requiring 
that all closed-loop eigenvalues lie within a specified sector. 

Figure 7a presents example step response envelopes and 
the response of the tip to a 4.8 cm position command input for 
500 Monte Carlo evaluations. The control time history 
corresponding to the mean response is given in Fig. 7b. The 
transfer function between hub torque and tip position is non- 
minimum phase; thus the step response exhibits an initial 
response in the wrong direction. The time response envelopes 
in Fig. 7 indicate the maximum acceptable non-minimum phase 
response. For 500 responses, the probability of violating the 
time response envelope is 0.184 with 95% confidence intervals 
(L,U) = (0.151,0.221). Individual responses characteristic of 
those evaluated by Monte Carlo analysis are given in Fig. 8. 
While the responses fill out the envelope, some of the 
individual responses within the envelope may not be acceptable 
in the face of real-world criteria governing rate of change of the 
response (Fig. 8c). This is a case where checking envelopes 
around the derivative of the response may be necessary. 
Similar analyses can be performed on control trajectories to 
make sure bandwidth and control effort limitations are not 
violated during the simulation. 

It is instructive from a design standpoint to plot stochastic 
robustness measures versus design parameters used to calculate 
feedback gains. Since there is only a single control in this 
example, the control weighting matrix R is a scalar and can be 
used as the design parameter. Two stochastic performance 
robustness measures are plotted versus R in Fig. 9 - the 
‘probability of violating the time response envelope and the 
probability of degree of instability. As control gains increase, 
the closed-loop roots are pushed farther into the left half plane, 
but they tend to migrate farther from their nominal values as 
well. At some value of control gain, there is a tradeoff 
between how far roots migrate and their location in the left half 
plane; thus a local minimum is apparent around R = 0.001 in 
the probability-of-degree-of-instability curves. While degree of 
instability improves for very small R, the control gains for this 
case are unrealistically large. For larger R (smaller control 
gains), the nominal closed-loop roots have real parts in the 
range of the values of Z used; thus the probability of degree of 
instability increases rapidly beyond R = 0.01. T he probability 
of the violating time response envelopes takes on the same 
shape as a function of R as the probability of degree of 
instability. In this example, the probability of instability is 
zero for all values of R checked, however, minima in the 
probability of instability versus design parameter curve can 
occur as well [11]. 

Figure 10 show's the stochastic root locus for the LQG 
system with estimator gains based on disturbance effect matrix 
L = G and disturbance and noise covariance matrices 


W = 


ro.oos 

01 

L o 

10J 


(23, 24 
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With estimator states added, the stochastic root locus changes 
in overall character from the full-state feedback case. Peaks are 
sharper, and the distribution along the real axis is less 
pronounced. In particular, note the eigenvalues associated 
with the largest peaks. In the full-state feedback system (Fig. 4) 
a broad distribution is associated with these eigenvalues, yet in 
Fig. 10, this pair of eigenvalues undergo very little variation 
from their nominal values! While the extent of the distribution 
into the left half plane is about the same as in Fig. 4, LQG 
system eigenvalues do migrate into the right half plane. The 
probability of instability estimate and 95% confidence intervals 

for 50,000 evaluations are P = 0.0771, and (L,U) = 0.0748, 
0.0795, representing a significant loss in the stability 
robustness characteristic of the LQR system. These results 
reflect the well known fact that guaranteed stability margins of 
an LQR system are lost when an estimator is added. 
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Loop Transfer Recovery (LQG/LTR) [12] is a common 
design technique by which the trade-off between estimator 
performance and stability robustness can be made 
systematically using a single design parameter. Figure 11 
illustrates the LTR mechanism for this example using the 
probability of instability as the robustness measure. Here 
plotting the probability of instability as a function of design 
parameter v (W = vW 0 ) shows that there is a value of v (v = 

2) that minimizes the probability of instability. The fact that 
such a minimum exists may not be apparent by simply 
examining the estimator eigenvalues. The kind of results 
presented in Fig. 11 offer design insight and show robustness 
characteristics not revealed by other robustness metrics. 

CONCLUSIONS 

Stochastic robustness offers a rigorous yet straightforward 
alternative to current metrics for control system robustness that 
is simple to compute and is unfettered by normally difficult 
problem statements, such as non-Gaussian statistics, products 
of parameter variations, and structured uncertainty. Principles 
behind stochastic robustness can be applied to scalar 
performance metrics as well as time responses, making it a 
good candidate for overall robustness analysis. Both 
performance and stability metrics resulting from stochastic 
robustness analysis can provide details relating robustness 
characteristics intrinsic to a given application and scalar control 
design parameters, making it a good candidate for optimization 
techniques as well. 
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Figure 5 Histogram and cumulative distribution for the single-link robot 
with uniform parameters, 50,000 evaluations. 
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Figure 6 Top view of a stochastic root locus with sector bounds for the 
single-link robot, 50,000 evaluations. 



Figure 8 Examples of individual tip responses 

a) Acceptable response within envelope. 

b) Response violates envelope. 

c) Response is within envelope, but criteria governing its derivative 
may be required. 



Figure 10 Stochastic root locus for the single-link robot with state 
estimation (LQG), 50,000 evaluations. 
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Figure 7 Time histories associated with tip position command of 4.8 cm. 

a) 500 Monte Carlo evaluations of the tip response. Envelopes arc 
defined by scalar performance criteria. Nominal response is indicated by 
the solid line. 

b) Nominal control input. 



R 

Figure 9 Stochastic performance robustness metrics vs. control weighting 
matrix R. 

a) Probability of violating time response envelopes. 

b) Probability of degree of instability, for values L along the real axis. 
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Figure 1 1 Probability of Instability vs. LQG/LTR design parameter v. 
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